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Harvesting  the  thermal  energy  stored  in  the  ambient  environment  provides  a  potential  sustainable 
energy  source.  Thermoelectric  power  generators  have  advantages  of  having  no  moving  parts,  being  dura¬ 
ble,  and  light-weighted.  These  unique  features  are  advantageous  for  many  applications  (i.e.,  carry-on 
medical  devices,  embedded  infrastructure  sensors,  aerospace,  transportation,  etc.).  To  ensure  the  efficient 
applications  of  thermoelectric  energy  harvesting  system,  the  behaviors  of  such  systems  need  to  be  fully 
understood.  Finite  element  simulations  provide  important  tools  for  such  purpose.  Although  modeling  the 
performance  of  thermoelectric  modules  has  been  conducted  by  many  researchers,  due  to  the  complexity 
in  solving  the  coupled  problem,  the  influences  of  the  effective  Seebeck  coefficient  and  carrier  density 
variations  on  the  performance  of  thermoelectric  system  are  generally  neglected.  This  results  in  an  over¬ 
estimation  of  the  power  generator  performance  under  strong-ionization  temperature  region.  This  paper 
presents  an  advanced  simulation  model  for  thermoelectric  elements  that  considers  the  effects  of  both 
factors.  The  mathematical  basis  of  this  model  is  firstly  presented.  Finite  element  simulations  are  then 
implemented  on  a  thermoelectric  power  generator  unit.  The  characteristics  of  the  thermoelectric  power 
generator  and  their  relationship  to  its  performance  are  discussed  under  different  working  temperature 
regions.  The  internal  physics  processes  of  the  TEM  harvester  are  analyzed  from  the  results  of  computa¬ 
tional  simulations.  The  new  model  presented  by  this  paper  is  more  conservative  than  the  traditional 
model  in  predicting  the  performance  of  TEM  power  generator  working  under  strong-ionization  region. 
By  accounting  for  the  mechanisms  that  are  typically  ignored  or  simplified  in  the  existing  models,  this 
new  model  provides  more  holistic  descriptions  on  the  thermoelectric  power  generator  behaviors  and 
therefore  potentially  will  further  improve  the  accuracy  in  the  computational  simulations. 

©  2014  Elsevier  Ltd.  Ah  rights  reserved. 


1.  Introduction 

Thermoelectric  module  (TEM)  converts  thermal  and  electrical 
energy  with  no  moving  parts.  It  is  light-weighted  and  durable. 
When  serving  as  a  power  source,  TEM  generates  electric  power 
through  utilizing  thermal  energy  and  is  usually  called  thermoelec¬ 
tric  generator  (TEG).  It  has  been  employed  as  power  sources  in 
wristwatches  [1],  satellites,  space  probes,  and  unmanned  remote 
facilities.  Thermoelectric  module  is  also  promising  to  provide  inno¬ 
vative  solutions  for  energy  applications  in  transportation  arena. 
These  include,  for  example,  harvesting  the  waste  heat  from  vehicle 
engine  and  exhaust  gas  [2],  utilizing  the  thermal  energy  contained 
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in  pavement  [3]  and  other  infrastructure  to  provide  electricity  sup¬ 
plies  for  sensors,  etc.  Design  and  optimization  of  TEM  performance 
relies  on  understanding  the  materials  properties  and  fundamental 
energy  conversion  mechanism. 

Analytic  models  have  been  developed  to  describe  the  behaviors 
and  performances  of  thermoelectric  devices.  To  obtain  analytic 
solutions,  simplifications  have  to  be  made,  i.e.,  the  Thomson  effects 
and  temperature  dependency  of  material  properties  are  generally 
neglected  [4-8].  There  are  also  efforts  in  refining  the  models  to 
include  the  Thomson  effects  [9]  and  the  temperature  dependency 
of  material  properties  [10-12].  However,  there  are  inherent  limita¬ 
tions  of  analytical  method  due  to  its  incapability  to  account  for 
complex  processes  and  materials  properties. 

Models  based  on  electrical  analogy  method  utilizes  the  mature 
knowledge  in  the  area  of  electric  circuits  for  thermal  field  analysis 
[13,14].  Finite  difference  software  SPICE  (Simulation  Program  with 
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Nomenclature 

TEM 

thermoelectric  module 

Greek  letters 

TEG 

thermoelectric  generator 

a 

electrical  conductivity  (S/m) 

j 

electric  current  per  unit  area  (A/m2) 

a 

Seebeck  coefficient  (V/K) 

J  drift 

drift  current  (A/m2) 

a 

effective  Seebeck  coefficient  (V/K) 

J  diffusion 

diffusion  current  (A/m2) 

Sot 

delta  Seebeck  coefficient  (V/K) 

q 

heat  transfer  rate  per  unit  area  (W/m2) 

K 

thermal  conductivity  at  zero  current  (W/(m  K)) 

n 

outward  unit  normal  vector 

Kc 

convective  heat  transfer  coefficient  (W/(m2  K)) 

1/ 

electrostatic  potential  (V) 

P 

electrical  resistivity  (Q  m) 

T 

temperature  (K) 

Pm 

material  density  (kg/m3) 

Ttop 

top  boundary  temperature  (K) 

P 

electrochemical  potential  (eV) 

T bottom 

bottom  boundary  temperature  (K) 

P 

chemical  potential  (eV) 

Cm 

heat  capacity  (J/(I<  kg) 

T 

Thomson  coefficient  (V/K) 

e 

unit  charge  (C) 

n 

free  electron  density  (m-3) 

Subscripts 

P 

free  hole  density  (m-3) 

drift 

drift 

k  o 

Boltzmann  constant  (J/K) 

diffusion 

diffusion 

h 

Planck  constant  (J  s) 

bottom 

bottom  edge 

Ef 

Fermi  level  (eV) 

top 

top  edge 

Ec 

bottom  edge  of  the  conduction  band  (eV) 

m 

material 

Ev 

top  edge  of  the  valence  band  (eV) 

c 

convective 

Ed 

energy  level  for  donors  (eV) 

C 

conduction  band 

Ea 

energy  level  for  acceptors  (eV) 

F 

Fermi 

Nc 

effective  density  of  states  in  conduction  band  (m-3) 

V 

valence  band 

Nv 

effective  density  of  states  in  the  valence  band  (m-3) 

D 

donors 

nd 

donor  impurity  concentration  (m-3) 

A 

acceptors 

Na 

acceptor  impurity  concentration  (m-3) 

n 

electron 

K 

electron  effective  mass  (kg) 

P 

hole 

m; 

hole  effective  mass  (kg) 

0 

constant 

m0 

electron  rest  mass  (kg) 

D 

carrier  diffusion  coefficient  (m2/s) 

Integrated  Circuit  Emphasis),  which  is  widely  used  in  the  area  of 
circuit  analysis,  has  been  used  for  modeling  the  thermoelectric 
harvesters  [15-19].  The  advantage  of  this  approach  is  that  the  ther¬ 
mal  field  and  electric  field  can  be  easily  coupled  in  the  same  sim¬ 
ulation  environment.  While  this  method  is  powerful  in  simulating 
complex  circuit  loads  connected  to  the  thermoelectric  power  gen¬ 
erator,  use  of  such  tools  requires  specialized  training.  This  can  pose 
significant  limits  on  its  applications.  The  electrical  abstraction  of 
the  thermoelectric  device  typically  focuses  on  the  influence  of 
the  macroscopic  performance  of  the  TEM,  such  as  the  output 
power,  and  temperature  difference  between  the  two  ends.  Param¬ 
eter  distributions  (especially  three  dimensional  distributions) 
inside  the  TEM  are  difficult  to  be  visualized  with  electrical  analogy 
method.  This  method  lacks  the  sensitivity  of  the  module  size  on 
the  TEM  performance,  which  makes  it  difficult  to  optimize  the 
shape  of  the  TEM.  In  order  to  increase  the  visualization  ability  of 
the  electrical  analogy  method,  a  three  dimensional  TCAD  (Synop- 
sys  Technology  Computer  Aided  Design)  implementations  have 
been  conducted  20-23].  However,  the  governing  equations  and 
the  performances  of  the  TEM  model  still  need  to  be  further  verified. 

With  the  advancement  in  its  multi-physics  simulation  capabil¬ 
ity,  the  finite  element  method  (FEM)  has  become  an  attractive 
method  to  simulate  the  TEM  performance  numerically  24-29]. 
Thomson  effects  and  temperature  dependent  properties  of  the 
TEM  materials  are  easily  coupled  into  the  governing  equations 
[30-32].  The  finite  element  method  features  the  advantages  of 
user-friendly  interface  for  model  construction  and  results  visuali¬ 
zation  as  well  as  accuracy  in  the  simulation  results  [33].  The 
multi-physics  simulation  capability  helps  to  couple  the  interac¬ 
tions  between  the  thermal,  electrical,  as  well  as  other  physics  fields 
such  as  the  mechanical  field  [34,35]. 


The  current  density  equation  used  by  these  previous  finite 
element  models  can  be  summarized  as  Eq.  (1 ).  Eq.  (1 )  does  not  con¬ 
sider  the  effects  of  chemical  potential  on  the  properties  of  thermo¬ 
electric  materials,  which  however  has  been  found  to  have  major 
impacts  on  the  thermoelectric  conversion.  The  chemical  potential 
can  be  induced  through  using  the  effective  Seebeck  coefficient, 
which  will  be  introduced  in  details  in  the  subsequent  context  of 
this  paper. 

J  =  -a .  (W  +  a  •  VT)  (1) 

where  the  vector  J  is  the  electric  current  per  unit  area;  o  is  the  elec¬ 
trical  conductivity;  V  is  the  electrostatic  potential;  a  is  Seebeck 
coefficient  and  T  is  the  temperature. 

Besides,  Eq.  (1)  is  only  applicable  for  homogeneous  conditions, 
where  carrier  density  is  uniform  in  the  semiconductor  materials 
that  form  the  legs  of  TEM.  However,  TEM  has  different  tempera¬ 
tures  at  two  ends  under  operation  conditions  in  order  to  form 
the  thermal  gradient.  Correspondingly,  carrier  (electrons  or  holes) 
densities  are  different  at  the  two  ends  under  non-strong-ionization 
regions  (including  low-temperature  weak-ionization  region  and 
high-temperature  intrinsic-ionization  region,  etc.)  due  to  the  dif¬ 
ferences  in  the  excitation  energies  under  different  temperatures. 
This  will  drive  carriers  to  diffuse  along  the  density  gradient.  Mahan 
[36  pointed  out  that  the  density  always  varies  as  long  as  there  is  a 
current  flowing  in  the  leg  and  only  by  introducing  variations  in  the 
particle  density  can  we  eliminate  the  violation  of  Poisson’s  equa¬ 
tion  under  homogeneous  conditions. 

This  paper  describes  an  improved  model  for  TEM  that  considers 
the  effects  of  chemical  potential  and  carrier  density  variations.  We 
firstly  derive  the  general  governing  equations  and  auxiliary  rela¬ 
tionships  that  describe  the  physical  processes  in  TEM,  using  n-type 
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semiconductor  material  as  an  example.  Chemical  potential  (i.e.,  the 
effective  Seebeck  coefficient)  is  discussed  under  homogeneous 
condition  and  carrier  density  variation  is  considered  under  inho¬ 
mogeneous  condition.  The  equation  system  is  solved  with  finite 
element  platform  COMSOL  for  a  typical  thermoelectric  generator 
unit.  The  characteristics  and  performances  of  the  module  are 
discussed. 


Plugging  Eq.  (6)  into  (4),  we  have 
J  — * .  W  -  *«  ■  VT  -  *(gt)  VT  -  »(*)  Vn  (7) 

2  A.  Effective  Seebeck  coefficient 


2.  Governing  equations 


The  equations  governing  the  multi-dimensional  temperature 
and  electrical  potential  distributions  in  thermoelectric  material 
are  the  energy  conservation  equation  and  charge  conservation 
equations  [37]  in  the  absence  of  an  applied  magnetic  field. 


pmCm^  =  W(KVT)  +  p]2 


TJ 


VT  +  (Va)T 


(2) 


VJ  =  0 

where 

=  -ff[v(^v)+aVl] 


J  =  -C 


v(  E 


-aVT 


(3) 

(4) 


q  =  aTJ  -  kVT  (5) 

Here  the  vector  J  is  the  electric  current  per  unit  area;  the  vector  q  is 
the  heat  transfer  rate  per  unit  area;  Cm  is  the  heat  capacity  with  unit 
of  J/(K  kg);  pm  is  the  material  density  with  unit  of  kg/m3;  T  is  the 
temperature;  k  is  the  thermal  conductivity  at  zero  current;  a=  1/ 
p  is  the  electrical  conductivity;  p  is  the  electrical  resistivity;  a  is 
the  Seebeck  coefficient;  ju  =  p  +  eV  is  the  electrochemical  potential; 
p  is  the  chemical  potential;  V  is  the  electrostatic  potential;  and  e  is 
the  unit  charge  of  charged  particles  giving  rise  to  the  electrical  cur¬ 
rent,  which  is  positive  for  hole  and  negative  for  electron.  All  mate¬ 
rial  properties  (a,  k ,  cr,  etc.)  are  functions  of  temperature.  The 
Seebeck  coefficient  can  also  be  dependent  on  other  factors,  such 
as  doping,  besides  temperature.  However,  they  are  neglected  for 
simplicity  in  this  paper. 

The  last  term  of  Eq.  (2)  includes  both  Peltier  effect  and  Thomson 
effect.  Peltier  effect  describes  heat  absorbed  by  the  junction  from 
the  environment  and  is  calculated  as  the  term  7J(Va)T.  At  the 
junction  interface,  where  two  different  materials  are  in  contact, 
temperature  is  the  same  across  the  junction  because  of  heat 
exchange.  Heat  absorption  is  due  to  the  differences  in  the  Seebeck 
coefficients  of  junction  materials.  Inside  the  TEM  element  where 
thermoelectric  material  remains  the  same,  Peltier  heat  equals  to 
zero. 

Thomson  effect  is  represented  by  the  term  TJ(|f)VT  =  tJ VT, 
where  t  is  the  Thomson  coefficient.  The  Thomson  effect  is  caused 
by  the  Seebeck  coefficient  gradient  which  results  from  the  temper¬ 
ature  gradient  along  the  thermoelectric  material. 

The  sign  in  front  of  the  term  TJ[(J£)VT  +  (Va)T]  of  Eq.  (2)  is 
dependent  on  the  sign  convention  of  Thomson  coefficient.  Under 
the  convention  that  the  Thomson  coefficient  t  is  considered  to  be 
positive  if  the  conductor  tends  to  become  cooler  (or  heat  is 
released)  when  the  current  flows  in  the  direction  of  the  positive 
temperature  gradient  [38],  the  sign  of  the  last  term  is  negative. 
Under  the  opposite  convention  where  t  >  0  implies  that  heat  is 
absorbed  when  current  flows  towards  regions  of  higher  tempera¬ 
ture  [6],  the  sign  in  front  of  the  last  term  is  positive.  Generally 
the  former  convention  is  used  by  most  researchers. 

The  chemical  potential  p  is  a  function  of  particle  density  and  the 
temperature  T  [39]. 


Similar  to  the  second  term  in  Eq.  (7),  the  third  term  also  corre¬ 
sponds  to  electrical  conductivity  a  and  temperature  gradient  VT. 
In  other  words,  the  factor  plays  a  role  on  the  Seebeck  coeffi¬ 
cient,  which  is  included  in  the  term  of  effective  Seebeck  coefficient 
as  in  Eq.  (8)  [36]. 

s=a+(<§0n  =  a+5a  (8) 

where  term  (§jjf)n  is  defined  as  the  delta  Seebeck  coefficient.  Then 
Eq.  (7)  becomes 

J  =  -a-VV  -aa-VT-afAf)  Vn  (9) 

\cOTl J  j 

The  effective  Seebeck  coefficient  contains  the  derivative  of  the 
chemical  potential  with  respect  to  the  temperature,  as  well  as 
the  delta  Seebeck  coefficient.  For  n-type  semiconductor  material, 
the  chemical  potential  p  is  equal  to  the  difference  between  the 
Fermi  energy  level  EF  and  the  bottom  edge  of  the  conduction  band 
Ec  [6]. 

P  =  Ef-Ec  (10) 

Assuming  that  Ec  is  not  dependent  on  temperature,  there  is 


r<9(EF  -  Ec)] 


edT 


n 


(ii) 


From  theories  of  semiconductor  physics,  the  following  expres¬ 
sion  of  Fermi  energy  level  is  obtained  following  the  Maxwell- 
Boltzmann  statistics  for  n-type  material  under  different  ionization 
regions  corresponding  to  different  temperature  levels  [40]. 


Low-temperature  weak-ionization  region : 


Ef 


(12) 


Room-temperature  strong-ionization  region : 

Ef  =  Ec  +  koT  In  (13) 


High-temperature  intrinsic-ionization  region : 


(14) 


where  k0  is  the  Boltzmann  constant  and  equals  to  1.38  x  10-23  J/K; 
Ed  is  the  energy  level  for  donors;  Ev  is  the  top  edge  of  the  valence 
band;  ND  is  the  donor  impurity  concentration;  Nv  is  the  effective 
density  of  states  in  the  valence  band;  Nc  is  the  effective  density  of 
states  in  conduction  band  and  can  be  described  as  follows: 


Nc 


n  (2nm;k0T^ 
h 3 


(15) 


(16) 


where  m*  and  m*  are  the  electron  effective  mass  and  hole  effective 
mass,  respectively;  h  is  Planck  constant. 
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Assuming  ED,  Ec ,  Ev,  ND,  m*  and  m*  are  independent  of  temper¬ 
ature  and  plugging  Eqs.  (12)-(16)  into  (11),  we  can  get  the  results 
of  delta  Seebeck  coefficient  as  follows: 


Weak-ionization  region  : 
Strong-ionization  region  : 
Intrinsic-ionization  region  : 


(17) 

(18) 
(19) 


From  Eqs.  (15)-(19),  for  the  weak  and  strong-ionization  regions,  the 
delta  Seebeck  coefficient  depends  on  the  donor  impurity  concentra¬ 
tion  (Nd),  the  electron  effective  mass  (m*),  and  the  temperature. 
The  relationships  can  be  referred  to  Fig.  1  using  the  strong-ioniza¬ 
tion  region  as  an  example. 

Fig.  1  shows  that  the  delta  Seebeck  coefficient  dot  is  positive  for 
n-type  material  under  strong-ionization  region  because  the  donor 
impurity  concentration  ND  is  less  than  the  effective  density  of 
states  in  conduction  band  Nc  and  e  is  negative  (referring  to  Eq. 
(18)).  Usually  the  Seebeck  coefficient  ot  is  negative  for  n-type  mate¬ 
rial.  Therefore,  the  signs  of  dot  and  ot  are  opposite  under  strong-ion¬ 
ization  region,  which  is  also  true  for  weak-ionization  region  of  both 
n-type  and  p-type  material,  as  well  as  strong-ionization  region  of 
p-type  material. 

For  the  intrinsic-ionization  region  in  the  case  of  n-type  material, 
dot  changes  as  the  ratio  between  the  hole  effective  mass  (m*)  and 
electron  effective  mass  (m*)  changes,  m*  and  m *  can  be  calculated 
according  to  the  corresponding  band  structure  of  the  thermoelec¬ 
tric  semiconductor  material.  When  m *  <  m*,  dot  has  an  opposite 
sign  of  ot.  When  m*  >  m*,  dot  has  the  same  sign  of  ot.  For  p-type 
material,  the  delta  Seebeck  coefficient  dot  is  the  same  as  n-type 
material.  Therefore,  when  m*  <  m*,  dot  has  the  same  sign  as  a; 
when  m*  >  m*,  dot  has  an  opposite  sign  of  ot. 

Besides,  the  Seebeck  coefficient  of  the  most  used  thermoelectric 
power  generator  material  is  of  several  hundreds  of  micro-volts  per 
Kelvin  (pV/K).  It  can  also  be  seen  from  the  figure  that  the  absolute 
value  of  dot  is  comparable  to  a,  especially  for  material  with  high 
effective  electron  mass  (m*)  under  low  doping  level.  This  results 
in  a  significant  influence  of  the  effective  Seebeck  coefficient,  as 


well  as  the  current  caused  by  temperature  gradient  (referring  to 
Eq.  (9)).  Therefore,  this  factor  should  not  be  neglected  as  Mahan 
[36  pointed  out.  The  effective  Seebeck  coefficient  is  considered 
in  the  model  implemented  in  this  paper. 

With  the  calculation  of  the  delta  Seebeck  coefficient  and  the 
effective  Seebeck  coefficient,  we  can  calculate  the  current  driven 
by  temperature  gradient  via  Eq.  (9).  Actually,  the  temperature  gra¬ 
dient  leads  to  a  voltage  potential  gradient  in  TEM  through  effective 
Seebeck  coefficient.  Therefore  current  driven  by  temperature  gra¬ 
dient  and  potential  gradient  can  be  included  as  drift  current  as 
follows: 

Idrift  =  •  VV  -  ca  •  VT  (20) 


2.2.  Carrier  density  variation 


The  last  term  cr(^)TVn  in  Eqs.  (7)  or  (9)  describes  diffusion  cur¬ 
rent  introduced  by  carrier  density  gradient  in  the  semiconductor 
thermoelectric  material.  Usually  in  semiconductor  physics,  the 
relationship  between  the  diffusion  current  density  and  carrier  den¬ 
sity  gradient  is  connected  through  the  carrier  diffusion  coefficient, 
following  Eq.  (21). 


J  diffusion 


=  eD  ■  Vn 


(21) 


Comparing  Eq.  (21)  with  term  0"(|^)TVn,  we  can  derive  the 


(22) 


expression  of  D: 

D  =  —  )  =  * 

e2\dn)T  e 2 

r9(£F-£c)i 

dn 

In  n-type  semiconductor  material,  the  density  of  carrier  elec¬ 
trons  is: 


(23> 

The  derivative  in  Eq.  (22)  is  shown  in  Eq.  (24)  under  the 
assumption  that  Nc  is  independent  of  n: 


d(Ep  —  Ec)  _  koT 

dn  ~  n 

Plugging  Eq.  (24)  into  (22),  we  have 


(24) 


-1.340E-04 

-5.625E-05 

2.150E-05 

9.925E-05 

1.770E-04 

2.548E-04 

3.325E-04 

4.103E-04 

4.880E-04 

5.658E-04 

6.435E-04 

7.213E-04 

7.990E-04 

8.768E-04 


Fig.  1.  The  relationship  between  the  delta  Seebeck  coefficient  and  the  donor  impurity  concentration  (ND),  the  relative  electron  effective  mass  (m*/m0)  and  the  temperature 
under  the  strong-ionization  region  for  n-type  thermoelectric  material,  where  m0  is  the  electron  rest  mass. 
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(T /'dfi\  _ak0T 
e2  \dn)T  ~  e2n 


(25) 


The  total  current  can  be  calculated  from  the  derived  the  expres¬ 
sions  of  drift  current  and  diffusion  current: 


J  J drift  T"  J diffusion 

J drift  =  -G  •  W  -  GOL  ■  VT  (26) 

J  diffusion  =  -eD  Vn 

Eq.  (26)  reveals  the  nature  of  current  flowing  inside  the  thermo¬ 
electric  material.  It  is  made  up  by  drift  current  and  diffusion  cur¬ 
rent.  The  drift  current  is  driven  by  the  electric  field  governed  by 
the  gradient  of  electric  potential  and  effective  Seebeck  potential. 
The  diffusion  current  is  caused  by  the  variations  in  the  carrier  con¬ 
centration  and  is  usually  neglected  by  most  current  models  in  the 
published  literature.  The  model  developed  in  this  paper  is  an 
improvement  over  existing  models  by  considering  the  effects  of 
both  drift  currents  and  diffusion  currents. 

IN  SUMMARY,  the  governing  equations  and  auxiliary  equations 
for  the  holistic  computational  model  of  TEM  are  listed  in  the 
following: 


Governing  equations : 


I  PmCm§=  V(KVT)  +  pi 2  -  Tj [(f) VT  +  ( Va)r] 
J  =  -<7  W  -  era  VT-eD  Vn 
<  J  =  —  <7  •  VV  —  era  •  VT  —  eD  •  Vp 
V-J  =  0 

k  q  =  aTJ  -  kVT 

(27) 


Auxiliary  equations  : 


r  n-type  material 

(2t im*nk0T)2 

Nc  =  2 - p - 

n=Nc  exp  (-^f) 


a  =  a  +  Sol 


p-type  material 

_  n  (2nm*pk0T)l 
V  ~  h3 

P  =  Nvex  p(^f) 


a  =  a  +  Sol 


(28) 


where  Nv  is  the  effective  density  of  states  in  valance  band.  EF  and  Sol 
can  be  referred  to  Table  1  corresponding  to  different  temperature 
levels. 


and  is  entirely  contributed  by  the  ionization  of  doped  impurity. 
Therefore  the  diffusion  current  (Eq.  (21))  is  zero.  This  scenario 
can  be  seen  as  the  so-called  homogeneous  case,  which  can  be 
applied  to  TEMs  used  in  the  room  temperature  range. 

If  the  semiconductor  thermoelectric  material  is  doped  in  gradi¬ 
ent,  or  working  under  weak-ionization  region  and  intrinsic-ioniza¬ 
tion  region  even  though  doped  evenly,  the  previous  phenomenon 
of  uniform  carrier  concentration  does  not  happen.  Carrier  concen¬ 
tration  increases  dramatically  as  temperature  increases.  In  such 
inhomogeneous  scenarios,  the  diffusion  current  arouses  and  is  dri¬ 
ven  by  the  carrier  concentration  gradient.  TEM  applications  under 
low-temperature  range  (for  example,  outer  space  applications)  and 
high-temperature  range  (for  example,  thermal  power  plant  appli¬ 
cations)  follow  these  phenomena. 

Fig.  2  gives  the  schematic  of  carrier  concentration  variations 
with  temperature  for  a  typical  n-type  thermoelectric  semiconduc¬ 
tor  material. 


3.  Holistic  modeling  of  the  TEM  performance,  results,  and 
discussions 

The  geometry  of  the  thermoelectric  generator  is  described  in 
the  following  Fig.  3.  The  semiconductor  legs  are  assumed  to  have 
dimensions  of  1.5  x  1.5  x  2  mm3.  Space  between  two  legs  is 
1  mm.  Metal  inter-connector  and  electrodes  are  assumed  to  have 
thickness  of  0.1  mm.  Both  the  open  circuit  and  close  circuit 
conditions  are  simulated.  The  resistor  is  assumed  to  be 
1.5  x  1.5  x  0.1  mm3  when  the  thermoelectric  generator  powers  a 
load. 

Properties  of  semiconductor  materials  are  described  as  temper¬ 
ature  dependent.  The  Seebeck  coefficient  of  n-type  material  is 
assumed  to  follow  Eq.  (30).  For  simplicity,  we  assume  that  the  See¬ 
beck  coefficient  of  p-type  material  has  the  opposite  sign  of  n-type, 
nevertheless,  more  complex  material  properties  can  be  easily 
incorporated.  The  electrical  conductivity  and  thermal  conductivity 
of  both  p-type  and  n-type  are  assumed  to  be  equal  and  follow  Eqs. 
(31)  and  (32),  respectively.  The  temperature  dependences  of  these 
parameters  can  also  be  referred  to  as  Fig.  4.  Other  material  param¬ 
eters  of  the  semiconductor  are  listed  in  Table  2. 

a  =  348.98  -  1.8557  + 0.0061 604T2  -  6.4013  x  10  673 

+  2.1847  x  1CT9T4  |.lV/K  (30) 


2.3.  Homogeneous  and  inhomogeneous  conditions 


a  =  70  x  (116.83  +  13561e-°0068336T)  S/cm  (31) 


For  strong-ionization  region,  if  we  plug  Eq.  (13)  into  (23),  we 
can  get  the  carrier  concentration  of  electrons  as 

n  =  ND  (29) 

When  the  semiconductor  thermoelectric  material  is  evenly 
doped  with  impurities,  the  carrier  concentration  inside  the  mate¬ 
rial  keeps  as  a  constant  everywhere  when  temperature  changes 


k  =  0.23912 


609.95 

T 


838.3  .....  ... 

~Y~  W/(mI<) 


(32) 


Material  properties  of  the  metal  inter-connector,  electrodes, 
and  wires  are  defined  as  constants  as  following:  the  Seebeck  coef¬ 
ficient  is  -3.5  pV/K;  the  electrical  resistivity  is  28.2  nQ  m;  the 
thermal  conductivity  is  205  W/(m  K);  the  heat  capacity  is  0.897  J/ 
(g  K);  and  the  density  is  2.7  g/cm3. 


Table  1 

Expressions  of  EF  and  Sot  under  different  temperature  levels  for  n-type  material. 


Low-temperature  weak-ionization  region 

Room-temperature  strong-ionization  region 

High-temperature  intrinsic-ionization  region 

n-Type 

Ef 

¥?a  +  ¥ln(:$() 

Ec  +  k0T ln 

&±fv  +  ^ln 

Sot 

¥[ta(&)-i] 

p-Type 

Ef 

Ev+Ea  k0T  ln  (  Na  \ 

2  IT111  \  2NV ) 

Ev-k0Tln(%) 

Sot 

¥  [1"  (»)-!] 

-$m@) 
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Fig.  2.  Rough  relationship  between  carrier  concentration  and  temperature  in 
semiconductor  thermoelectric  material  impurity. 


Finite  element  tool  COMSOL@  is  utilized  to  develop  the  model 
and  solve  the  PDEs  numerically.  Coefficient  form  PDE  modules 
are  used  to  describe  the  governing  equations  and  boundary  condi¬ 
tions.  Auxiliary  equations  are  imported  as  functions.  The  shape  of 
elements  are  triangle  on  each  surface  and  thus  in  tetrahedron 
shape  in  three  dimensions.  Each  node  has  one  degree  of  freedom 
for  temperature,  one  degree  of  freedom  for  voltage  potential  and 
three  degrees  of  freedom  for  current.  A  grid  sensitivity  analysis  is 
performed  on  the  sensitivity  of  results  on  the  size  of  FEM  mesh. 
The  results  of  four  different  FEM  mesh  options  are  compared. 
The  simulation  case  is  a  TEM  working  under  open  circuit  condition. 
The  overall  output  voltages  are  monitored  when  the  bottom  sur¬ 
face  temperature  of  the  TEM  is  kept  at  340  K  and  the  upper  surface 
temperature  ranges  from  300  K  to  330  K.  The  results  are  summa¬ 
rized  in  Table  3.  Physics-controlled  fine  mesh  is  selected  in  the 
subsequent  simulations  based  on  a  comprehensive  consideration 
of  accuracy  and  computational  time. 


Metal 

Inter-connector 


(a) 


Fig.  3.  Geometry  of  the  thermoelectric  generator  (a)  open  circuit  (b)  with  a  load  resistor. 


Fig.  4.  Temperature  dependence  of  material  properties  (a)  Seebeck  coefficient  of  n-type  material  (b)  electrical  conductivity  and  thermal  conductivity. 


Table  2 

Material  properties  of  semiconductor  thermoelectric  material. 


Parameter 

Symbol 

Value 

Unit 

Parameter 

Symbol 

Value 

Unit 

Density 

Pm 

11 

g/cm3 

Heat  capacity 

Cm 

0.1625 

J/(gK) 

Electron  effective  mass 

mn 

0.4mo 

kg 

Hole  effective  mass 

mp 

1.1  m0 

kg 

Bottom  edge  of  conduction  band 

Ec 

0.67 

eV 

Top  edge  of  valence  band 

Ey 

0 

eV 

Energy  level  for  donors 

Ed 

0.66 

eV 

Energy  level  for  acceptors 

Ea 

0.1 

eV 

Donor  impurity  concentration 

nd 

5  x  1018 

cm-3 

Acceptor  impurity  concentration 

Na 

5  x  1018 

cm-3 
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Table  3 

Grid  sensitivity  test  when  the  TEM  is  in  open  circuit. 


Open  circuit  overall  output  voltage 

Finer  mesh 

Fine  mesh 

Normal  mesh 

Coarse  mesh 

300  K 

4.96  mV 

4.96  mV 

4.95  mV 

4.95  mV 

Upper  surface  temperature 

310  K 

2.84  mV 

2.84  mV 

2.84  mV 

2.83  mV 

320  K 

1.29  mV 

1.29  mV 

1.29  mV 

1.29  mV 

330  K 

333.429  pV 

333.722  pV 

334.396  pV 

335.051  pV 

Approximate  computation  time 

110s 

74  s 

62  s 

51  s 

3A.  Behaviors  of  TEM  unit  under  strong-ionization  region 

Most  thermoelectric  power  generator  are  designed  to  work 
under  strong  ionization  region  which  corresponds  to  the  highest 
ZT  value.  As  methodologies  for  analyzing  TEM  performance  under 
weak,  strong,  and  intrinsic  ionization  regions  are  similar,  the 
strong-ionization  case  is  firstly  presented  as  an  example.  Compar¬ 
isons  between  different  ionization  regions  are  then  presented. 

For  the  strong-ionization  case,  the  bottom  boundaries  of  both  p- 
leg  and  n-leg  are  set  at  340  K.  Meanwhile  the  left  side  boundary  of 
the  bottom  left  metal  electrode  is  set  at  potential  of  0  V,  referring 
to  Fig.  3.  The  side  boundaries  of  all  solid  domains  are  assumed  to 
have  convective  heat  exchange  with  the  ambient,  following  Eq. 
(33). 

-n  •  (kVT)  =  Kc{Tbottom  -  5  -  T)  (33) 

where  n  is  the  outward  unit  normal  vector  seen  from  the  inside  of  a 
certain  solid  domain;  Tbottom  is  the  bottom  boundary  temperature 
and  T  is  the  temperature  at  the  solid  boundary.  Eq.  (33)  means 
the  ambient  temperature  is  assumed  to  be  Thottom  -  5  and  the  con¬ 
vective  heat  transfer  coefficient  kc  between  the  ambient  and  the 
module  is  assumed  to  be  10  W/(m2  I<).  When  the  environment  tem¬ 
perature  is  higher  than  the  boundary  temperature,  heat  is  absorbed 
into  the  module.  When  the  environment  temperature  is  lower  than 
the  boundary  temperature,  heat  flux  flows  outwards  as  seen  from 
the  solid  domain. 

3AA.  Open  circuit 

Fig.  5  shows  the  directions  of  current  driven  by  Seebeck  effect 
caused  e-field  and  delta  Seebeck  effect  caused  e-field  when  the 
thermoelectric  generator  works  under  10  K  temperature  gradient 
in  the  strong-ionization  region  (temperatures  at  the  top  surface 
and  bottom  surface  of  the  TEM  metal  electrode  are  330  K  and 
340  K  respectively). 

Under  this  temperature  range,  the  Seebeck  coefficients  of  the  p- 
type  material  and  n-type  material  are  positive  and  negative  respec¬ 
tively.  Seebeck  coefficient  is  defined  as  the  following  equation, 


a  =  - 


AV seebeck 

AT 


(34) 


which  means  that  the  temperature  gradient  is  in  the  opposite  direc¬ 
tion  of  the  Seebeck  caused  e-field  for  positive  Seebeck  coefficient 
and  it  is  reverse  for  negative  Seebeck  coefficient.  As  the  tempera¬ 
ture  gradient  goes  downwards  for  both  p-leg  and  n-leg,  e-field 
caused  by  Seebeck  effect  points  upwards  in  the  p-leg  and  down¬ 
wards  in  the  n-leg,  as  is  shown  in  Fig.  5a.  It  has  been  pointed  out 
in  the  early  context  that  the  delta  Seebeck  coefficient  has  the  oppo¬ 
site  sign  of  Seebeck  coefficient  in  strong-ionization  region,  for  both 
p-type  material  and  n-type  material.  Whence  the  e-field  induced  by 
the  delta  Seebeck  coefficient  is  reverse  to  Seebeck  e-field,  as  shown 
in  Fig.  5b.  Under  the  material  property  definition  of  the  TEM  ana¬ 
lyzed  here,  the  absolute  value  of  delta  Seebeck  coefficient  is  much 
smaller  than  the  Seebeck  coefficient,  so  the  overall  effective  Seebeck 
e-field  is  still  in  the  same  direction  of  the  Seebeck  e-field,  as  shown 
in  Fig.  5c. 

Fig.  6  illustrates  the  directions  of  current  driven  by  voltage- 
potential-induced  e-field  and  carrier  density  gradient.  The  total 
current  is  the  superposition  of  those  current  components,  based 
on  Eq.  (9).  For  TEM  with  evenly  doped  impurities  and  working 
under  strong-ionization  condition,  the  carrier  density  is  constant 
everywhere  according  to  Eq.  (29).  This  results  in  a  zero  net  diffu¬ 
sion  current  density  as  shown  in  Fig.  6b.  Driven  by  the  effective 
Seebeck  e-field,  current  flows  from  the  p-leg  to  the  n-leg,  as  shown 
in  Fig.  6c.  A  positive  voltage  potential  forms  at  the  right  side  of  the 
bottom  right  metal  electrode,  which  cause  a  reverse  voltage  poten¬ 
tial  through  the  TEM,  as  depicted  in  Fig.  6a. 

Fig.  7  presents  the  distributions  of  temperature,  current  density, 
and  voltage  inside  the  thermoelectric  module.  Fig.  7a  shows  that 
the  temperature  changes  from  330  K  to  340  K  gradually  from  top 
surface  to  bottom  surface.  Eq.  (3)  used  in  the  simulation  model 
(i.e.,  V  •  J  =  0),  implies  that  there  is  no  current  source  or  sink  inside 
the  module.  It  forces  all  carriers  in  the  semiconductor  legs  to  flow 
through  the  metal  inner-connector.  Fig.  7b  shows  that  the  current 
magnitude  is  much  higher  in  the  top  metal  inner-connector.  Mean¬ 
while  the  current  amplitude  remains  relatively  stable  specifically 


Fig.  5.  Directions  of  current  driven  by  (a)  Seebeck  e-field  (b)  delta  Seebeck  e-field  and  (c)  effective  Seebeck  e-field  when  the  thermoelectric  generator  is  working  under  10  K 
temperature  difference  and  in  open  circuit  (arrows  only  illustrate  the  direction  of  current,  not  their  magnitudes). 
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Fig.  6.  Directions  of  current  driven  by  (a)  voltage  potential  e-field,  (b)  carrier  density  gradient  and  (c)  the  total  current  direction,  when  the  thermoelectric  generator  is 
working  under  10  K  temperature  difference  and  in  open  circuit  (arrows  only  illustrate  the  direction  of  current,  not  magnitudes). 


Fig.  7.  Distribution  of  (a)  temperature  (b)  current  density  and  (c)  voltage  when  the  thermoelectric  generator  is  working  under  10  K  temperature  difference  and  in  open 
circuit. 


for  the  p-leg  and  n-leg.  Fig.  7c  proves  that  the  output  voltage  is 
positive  when  temperature  gradient  is  going  upwards. 

When  the  temperature  of  top  boundary  of  the  metal  inter-con¬ 
nector  is  swept  from  300  K  to  340  K,  while  the  temperature  of  the 
bottom  boundary  of  the  metal  electrodes  remain  at  340  K,  the  one 
dimensional  temperature,  current  density  amplitude,  Thomson 
heat  and  voltage  distribution  along  the  middle  line  of  p-leg  and 
n-leg  (also  including  the  bottom  metal  electrodes  and  top  metal 
inner-connector)  are  shown  in  Fig.  8  respectively.  The  part  from 
0  to  0.1  mm  along  the  horizontal  axis  corresponds  to  the  bottom 
metal  electrodes.  The  part  from  0.1  mm  to  2.1  mm  corresponds 
to  the  p-leg  or  n-leg.  The  part  from  2.1  mm  to  2.2  mm  corresponds 
to  the  top  metal  inner-connector. 

From  the  material  properties  definitions  of  the  TEM  used  in  this 
simulation,  m*  >  m*.  So  NV>NC  according  to  Eq.  (29).  Considering 
Na  =  Nd  here,  so  we  know  that  the  absolute  value  of  delta  Seebeck 
coefficient  of  the  p-leg  is  larger  than  n-leg,  based  on  parameters 
listed  in  Table  1.  Therefore,  the  absolute  value  of  effective  Seebeck 
coefficient  of  the  p-leg  is  smaller  than  n-leg,  which  results  in  that 
the  current  density  amplitude  of  the  p-leg  is  smaller  than  that  of 
n-leg,  as  shown  in  Fig.  8b. 

The  Thomson  heat  by  the  semiconductor  material  is 
-TJ[(|f)VT  +  (Va)T],  according  to  Eq.  (2).  In  the  temperature  range 
corresponding  to  strong-ionization  region,  is  negative  for  p-leg 
and  positive  for  n-leg.  Because  J  and  VThave  the  opposite  direction 
in  p-leg  and  the  same  direction  in  n-leg,  the  Thomson  heat  term  is 
negative  for  both  p-leg  and  n-leg,  which  means  that  heat  is 
extracted  from  the  semiconductor  material.  The  higher  the  current 
density,  the  more  heat  is  extracted.  Fig.  8c  verifies  the  prediction, 
where  Thomson  heat  is  negative  for  both  n-leg  and  p-leg  and 
n-leg  absorbs  more  Thomson  heat.  It  can  also  be  noticed  that  there 
is  a  jump  in  the  absorbed  heat  around  the  junctions  between  semi¬ 
conductor  legs  and  the  metal  electrodes.  This  is  possibly  due  to  the 
Peltier  heat  change  in  addition  to  the  Thomson  heat. 


Even  though  n-leg  has  a  higher  current  density  than  p-leg 
(which  leads  to  a  higher  heat  loss  of  the  n-leg),  the  temperature 
distributions  of  p-leg  and  n-leg  are  approximately  the  same,  seen 
from  Fig.  8a.  This  is  because  more  Joule  heat  is  produced  inside 
the  n-leg,  which  compensates  its  relatively  higher  heat  loss. 

Fig.  8d  shows  the  open  circuit  voltage  potential  distributions 
under  different  temperature  differences,  where  the  dash  line  val¬ 
ues  are  the  overall  output  voltage  of  the  module.  It  implies  that 
the  relationship  between  the  output  voltage  and  the  temperature 
difference  across  the  TEM  is  nonlinear.  The  larger  the  temperature 
difference,  the  higher  the  slope  of  the  output  voltage  versus 
temperature  difference  curve,  i.e.,  the  higher  the  effective  Seebeck 
coefficient. 

3.1.2.  TEM  performance  with  a  resistive  load 

Simulations  are  also  conducted  when  the  TEM  is  powering  an 
electric  load  (i.e.,  a  resistor),  with  set  up  as  shown  in  Fig.  3b. 
The  resistivity  of  the  load  resistor  is  swept  from  10-8  to 
10-3Qm,  approximately  ranging  from  short  circuit  to  open  cir¬ 
cuit  conditions.  The  output  current  and  output  voltage  are  shown 
in  Fig.  9a  as  functions  of  load  resistance.  The  output  power  is 
obtained  from  the  output  current  and  voltage.  The  output  power 
efficiency  is  determined  by  further  dividing  the  absorbed  heat 
from  the  bottom  boundary  of  the  TEM.  The  variations  of  absorbed 
heat  and  the  output  power  efficiency  with  load  resistance  are 
shown  in  Fig.  9b. 

Fig.  9  implies  that  output  current  decreases  and  output  voltage 
increases  as  load  resistivity  goes  up.  The  absorbed  heat,  however, 
does  not  change  much  as  load  resistance  changes.  The  maximum 
output  power  efficiency  occurs  at  resistivity  of  round  1  pQ  m  (or 
output  resistance  of  0.01  Q),  which  corresponds  to  the  point  that 
the  load  resistance  is  equal  to  the  inner  resistance  of  the  TEM. 
Results  of  the  simulations  here  indicate  that  the  larger  the  temper¬ 
ature  difference  across  the  TEM  element,  the  larger  the  output 
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Fig.  8.  One  dimensional  distributions  of  (a)  temperature  (b)  current  density  amplitude  (c)  absorbed  Thomson  heat  and  (d)  voltage,  along  the  middle  line  of  p-leg  and  n-leg. 
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Fig.  9.  (a)  Functions  of  output  current  and  output  voltage  with  respect  to  different  load  resistivity  (solid  lines  represent  output  current.  Dash  lines  represent  output  voltage. 
Red  arrows  means  when  the  load  resistivity  goes  down  to  10-8  Q  m,  the  circuit  is  similar  to  a  short  circuit  condition.  When  the  load  resistivity  goes  up  to  10-3  Q  m,  the  circuit 
approaches  to  the  open  circuit  condition);  (b)  monitored  absorbed  heat  from  the  bottom  boundary  of  the  TEM  and  the  output  power  efficiency  as  functions  of  load  resistivity 
(solid  lines  correspond  to  power  efficiency.  Dash  lines  correspond  to  absorbed  heat).  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  this  article.) 


power  and  power  efficiency,  which  are  consistent  with  general 
observations  in  TEM. 

3.2.  Behaviors  of  TEM  unit  working  under  different  ionization  regions 
predicted  by  the  new  model  versus  traditional  model 

Simulations  are  also  conducted  when  TEM  unit  works  under 
weak-ionization  region,  strong-ionization  region,  and  intrinsic- 


ionization  region.  For  all  these  conditions,  temperature  difference 
of  10  K  is  assumed  across  the  TEM  element.  The  reference 
temperature  at  the  bottom  conductor  are  140  K,  340  K,  and 
1000  K  respectively.  The  same  boundary  conductive  heat  exchange 
is  assumed  (Eq.  (33)).  Computational  simulations  based  on 
traditional  model  following  Eq.  (1 )  are  also  conducted  as  references 
under  corresponding  temperature  conditions  (which  are  plotted  in 
dash  lines). 
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Fig.  10.  Open  circuit  thermoelectric  power  generator  performances  comparison  between  different  ionization  regions  when  the  TEM  is  working  under  10  K  temperature 
difference,  with  traditional  model  results  as  references,  including  one  dimensional  distributions  of  (a)  current  density  amplitude,  (b)  temperature  (portion  of  the  curves  are 
zoomed  in),  (c)  absorbed  Thomson  heat,  (d)  voltage  potential. 


Fig.  10a  shows  that  when  working  under  weak-ionization 
region,  the  current  amplitude  is  much  bigger  in  the  bottom  left 
metal  electrode  than  anywhere  else,  while  traditional  model  does 
not  predict  this  behavior.  For  weak  and  intrinsic  ionization  regions, 
the  new  model  predicted  higher  current  densities  than  the  tradi¬ 
tional  model  by  6%  and  195%.  This  is  due  to  the  big  diffusion  cur¬ 
rent  under  inhomogeneous  conditions  that  is  captured  by  this 
model.  In  strong-ionization  region,  the  new  model  predicts  a  32% 
smaller  current  density  than  the  traditional  model,  because  the 
delta  Seebeck  coefficient  is  opposite  to  that  of  the  Seebeck  effect 
under  strong-ionization  region.  It  can  also  be  seen  from  the  figure 
that  under  the  same  temperature  difference,  the  higher  the  work¬ 
ing  temperature  of  the  module,  the  lower  the  current  density 
amplitude.  This  result  is  possibly  due  to  the  fact  that  the  electrical 
conductivity  of  a  material  decreases  with  increasing  temperature. 

Fig.  10b  shows  the  change  of  the  relative  temperature  differ¬ 
ence  along  the  p-leg  and  n-leg.  The  relative  temperature  difference 
is  calculated  as  the  absolute  value  of  the  temperature  at  a  certain 
point  minus  the  bottom  boundary  temperature.  The  figure  indi¬ 
cates  that  the  relative  temperature  differences  are  approximately 
the  same  under  various  conditions.  The  only  remarkable  phenom¬ 
ena  appears  for  the  case  of  the  relative  temperature  difference  at 
p-leg  under  weak-ionization  condition,  which  is  lower  than  the 
conditions.  It  means  the  temperature  inside  the  p-leg  is  relative 
lower  than  under  other  conditions.  Referring  to  the  material  prop¬ 
erty  definitions  of  the  TEM,  the  term  is  negative  for  n-leg  and 
positive  for  p-leg,  which  is  reverse  to  the  scenario  under  strong 


and  intrinsic  ionization  regions.  Consequently,  heat  is  absorbed 
from  the  ambient  inwards  to  the  semiconductor  leg.  This  explains 
the  relative  temperature  difference  in  the  p-leg  under  weak-ioni¬ 
zation  region  is  lower  than  other  cases,  which  can  also  be  proved 
by  Fig.  10c  as  the  Thomson  heat  for  the  p-leg  under  weak-ioniza¬ 
tion  range  is  significantly  larger  than  those  under  the  other  cases. 
Such  effects,  however,  is  not  captured  by  the  traditional  simulation 
model. 

Fig.  lOd  plots  the  voltage  distribution  in  the  p-leg  and  n-leg  of 
the  TEM  unit  under  different  working  conditions.  The  results  imply 
that  the  current  direction  in  the  p-leg  under  weak-ionization  tem¬ 
perature,  which  goes  from  bottom  to  top,  is  opposite  to  those  in  all 
the  other  cases  where  the  current  goes  from  top  to  bottom.  Due  to 
the  large  magnitude  produced  voltage  in  the  p-leg,  the  resultant 
overall  output  voltage  of  the  module  is  negative.  Another  impor¬ 
tant  observation  from  this  figure  is  that  the  traditional  simulation 
model  implies  under  the  same  temperature  difference  of  10  K,  the 
higher  the  working  temperature  of  the  module,  the  lower  the  out¬ 
put  voltage.  However,  the  new  model  predicts  a  reverse  trend.  The 
results  of  the  new  model  indicates  that  as  the  working  tempera¬ 
ture  increases,  the  absolute  value  of  the  output  voltage  will  first 
decrease  and  then  increase.  For  temperatures  in  the  strong-ioniza¬ 
tion  region,  the  voltage  output  predicted  by  the  new  model  is  39% 
lower  than  that  predicted  by  the  traditional  model. 

Results  for  the  case  with  a  load  resistor  are  shown  in  Fig.  11a 
and  b.  The  absolute  value  of  output  current  decreases  and  output 
voltage  increases  when  load  resistivity  increases.  The  resistivity 
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Fig.  11.  Close  circuit  TEM  performance  comparison  between  different  ionization  regions  when  TEM  is  working  under  10  K  temperature  difference,  with  traditional  model 
results  as  references  (a)  output  voltage  and  current  as  functions  of  load  resistivity;  (b)  absorbed  heat  and  output  power  efficiency  as  functions  of  load  resistivity. 


corresponding  to  the  maximum  output  power  efficiency  increases 
as  the  temperature  of  the  TEM  increases.  It  is  about  0.01  pQ  m, 
1  pQ  m  and  10  pQ  m  for  the  weak,  strong,  and  intrinsic  ionization 
regions,  respectively.  This  might  be  related  to  the  increases  in  the 
electrical  resistance  of  TEM  materials  with  temperature.  For  the 
strong-ionization  region,  the  maximum  output  power  efficiency 
predicted  by  the  new  model  predicts  is  65%  lower  than  the  tradi¬ 
tional  model.  Meanwhile  the  output  current  and  voltage  corre¬ 
sponding  to  the  maximum  power  efficiency  working  point 
predicted  by  the  new  model  are  lower  than  the  traditional  model 
prediction.  This  is  because  that  the  delta  Seebeck  coefficient  is 
opposite  to  that  of  the  Seebeck  effect  when  working  under  temper¬ 
ature  in  the  strong-ionization  region. 

4.  Conclusion 

This  paper  describes  the  development  of  the  governing  equa¬ 
tions  and  auxiliary  relationships  to  implement  holistic  simulations 
of  TEM  behaviors.  The  thermoelectric  processes  of  Seebck,  Peltier, 
and  Thomson  effects  and  Joule  heat  are  included  in  the  model. 
Particularly,  this  model  considers  the  influence  of  the  effective 
Seebeck  coefficient  and  the  carrier  concentration  variations.  Differ¬ 
ent  working  temperature  regions  corresponding  to  transitions  in 
carrier  concentrations  are  also  considered.  Finite  element  simula¬ 
tions  are  conducted  on  a  thermoelectric  generator  unit  under  the 
open  and  close  circuit  conditions. 

For  the  open  circuit  scenario,  the  physics  phenomena  inside  the 
thermoelectric  module  during  thermoelectric  generation  process 
are  discussed.  In  order  to  analyze  the  movements  of  the  carriers 
inside  the  thermoelectric  power  generator,  current  driven  by  the 
thermal  gradient  is  decomposed  into  components  driven  by 
Seebeck  e-field,  delta  Seebeck  e-filed,  effective  Seebeck  e-field, 
voltage  potential  e-field  and  carrier  density  gradient,  respectively. 
Comparison  indicates  that  the  carrier  movements  are  mainly 
decided  by  the  Seebeck  effect  induced  electric  field,  even  though 
they  are  impaired  by  the  effective  Seebeck  effect  and  the  output 
voltage  induced  electric  field.  The  carrier  movements  because  of 
carrier  density  variation  under  strong-ionization  region  does  not 
play  a  major  role  when  the  thermoelectric  legs  are  evenly  doped. 
However,  they  have  significant  contribution  to  the  overall  current 
when  working  under  weak  and  intrinsic  ionization  regions.  This 
proves  that  the  density  variations  of  carriers  on  the  performance 
of  TEM  generator  cannot  be  neglected.  The  overall  temperature, 
current  density,  Thomson  heat  and  voltage  potential  distributions 
are  described.  The  results  indicate  that  the  relationship  between 
the  output  voltage  and  the  temperature  difference  across  the 
thermoelectric  power  generator  is  non-linear,  partly  due  to  the 


effects  of  temperature  dependency  of  thermoelectric  material 
properties. 

For  the  close  circuit  scenario  with  a  load  resistor,  simulation 
results  on  the  power  performances  of  the  thermoelectric  power 
generator  are  discussed.  The  performance  of  the  TEM  at  different 
temperature  ranges  (corresponding  to  different  ionization  regions) 
are  compared.  Results  of  traditional  model  are  introduced  as  the 
reference  base.  The  output  voltage,  current,  absorbed  heat  and 
power  efficiency  are  described.  The  power  generator  works  at  its 
best  efficiency  when  the  load  resistance  equals  to  the  inner  resis¬ 
tance.  The  results  show  that  the  effective  Seebeck  coefficient  and 
the  carrier  concentration  variations  have  significant  influences  on 
the  performance  of  the  TEM  elements.  When  the  thermoelectric 
power  generator  is  working  under  strong-ionization  region,  the 
output  voltage,  current  and  output  power  efficiency  predicted  by 
the  new  model  is  lower  than  those  predicted  by  the  traditional 
model.  As  the  traditional  model  has  been  proved  to  overestimates 
the  performance  of  the  power  generator  under  strong-ionization 
region  [29],  it  implies  that  the  new  model  has  the  potential  to 
describe  the  TEM  behaviors  more  precisely.  The  results  also  show 
that  ignoring  the  effects  of  effective  Seebeck  coefficient  (delta  See¬ 
beck  coefficient)  and  carrier  concentration  variations  as  commonly 
done  in  the  traditional  simulations  models  fails  to  capture  many 
unique  behaviors  of  TEM  module. 

This  holistic  simulation  model  provides  important  tool  to  assist 
the  design  optimization  and  performance  prediction  of  innovative 
TEM  applications. 
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